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A numerical method for Mean Field Games on networks 


Simone Cacace* Fabio Camilli' and Claudio Mar chi* 


Abstract 

We propose a numerical method for stationary Mean Field Games defined on a 
network. In this framework a correct approximation of the transition conditions at 
the vertices plays a crucial role. We prove existence, uniqueness and convergence of 
the scheme and we also propose a least squares method for the solution of the discrete 
system. Numerical experiments are carried out. 
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1 Introduction 

The Mean Field Game (MFG in short) theory has been introduced in flGl fT8j j to describe 
the limit behavior of differential games when the number of agents becomes very large. 
Models based on this theory can be used to investigate crowd dynamics, consensus forma¬ 
tion and various economical and social problems (as growth theory, environmental policy 
and formation of volatility in financial markets) in which the strategy of the single agent 
determines a collective behavior of the population (see mmm)- 

From a mathematical point of view, MFG theory leads to the study of a coupled 
system of two differential equations: a Hamilton-Jacobi-Bellman equation and a Fokker- 
Planck equation, describing respectively the optimal behavior of each single agent and the 
evolution of the whole population. There is a rapidly increasing literature concerning both 
the theoretical aspects and the applications of MFG (see the review paper [13]). 

A crucial point to extend the theory of MFG systems to networks is to find the 
appropriate transition conditions at the vertices in order to obtain a well posed mathe¬ 
matical problem, coherent with the applications. In j8j, it was considered a MFG system 
with quadratic Hamiltonian which, by an appropriate change of variable, can be trans¬ 
formed into a linear system of differential equations coupled only via the initial datum. 
A general class of stationary MFG systems on networks is considered in [7], where it is 
proved existence and uniqueness of classical solutions to the problem 

— vd 2 u + H(x, du) + A = V[m], 

< ud 2 m + d(m H p (x,du)) = 0, x E T (1.1) 

f r m(x)dx = 1, / r u(x)dx = 0. 
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Here the network T = (V, £) is a finite collection of points V := in M” indexed by I, 

connected by continuous, non self-intersecting arcs £ := {ej}j £ j indexed by J. Moreover, 
v = (v.j)j£j are strictly positive numbers, the Hamiltonian If is a collection 
where Hj are continuous, convex Hamiltonians defined on the arcs ej, and H p (x, q ) is the 
differential of p H y H(x,p) at p = q. Let us stress that H may be discontinuous at the 
vertices. 


The equations in (1.1), which have the same interpretation as in the classical MFG theory, 
are defined in terms of the coordinate parametrizing the arc e.j, j E J, and have to be 
complemented with appropriate conditions at the vertices. At each internal vertex v t we 
consider the transition conditions 


Y. Vjdjii(vi ) = 0, 


felnci 


Y Wjdjm(vi) + Hj !P (vi,dju)mj(vi)] = 0, 

jelnci 

Uj{vi) = Uk(vi), rrijivi) = m k (vi), j, k E Inc* 


( 1 . 2 ) 


where Inc* denotes the set of the edges incident the vertex Vi and Hj tP is the derivative 
with respect to p of the Hamiltonian defined on the edge ej. We mention, see [7] for more 


details, that the first condition in (1.2) is the classical Kirchhoff condition and it prescribes 


the probability that an agent reaching the vertex V{ enters in the incident edge ej, j E Inc*; 


the second condition in (1.2) guarantees the mass conservation at v\ (the sum of the fluxes 
at Vi is null); the third condition is the continuity of u and m at v%. We remark that 
(1.2) are natural conditions for 2 nd order problems on networks. In fact the 


domain of the Laplace operator on a network is given by continuous functions 
on T which are H 2 on the edges and which satisfy the Kirchhoff condition at 


the vertices 


Moreover, the transition conditions are a crucial ingredient 


for the validity of the maximum principle on networks. 


In this paper we consider the numerical approximation of the problem (1.1)-(1.2) 


following the approach in m, where a finite difference approximation of the MFG system 
is studied (see also ia naira for different approaches). Inside the edges we follow the same 


approach of |2j and we discretize the differential equations in (1.1) by finite differences 


The guideline to find the correct approximation of the transition conditions in (1.2) is 


to reproduce at a discrete level some fundamental identities which are obtained in the 
continuous setting by the weak formulation of the problem (see f.e. (3.22)). For this 


reason the discrete Hamiltonian defined by a monotone approximation of the Hamilton- 
Jacobi-Bellman equation is also used in the discretization of the Fokker-Planck equation 
and of the corresponding transition condition. By means of the previous identities we 
prove the well-posedness of the discrete problem and the convergence to the solution of 


the system (1.1). 


While there is a large literature about the approximation of hyperbolic 
problems on networks (see for example [5j, [10]), as far as we know, numerical 
schemes for second order differential equations on networks with Kirchhoff 
conditions have been only considered in the linear case (see [ 19 ], [ 20 ] ). Hence 
the part concerning the approximation of the Hamilton-Jacobi-Bellman equation on the 
network is new and of independent interest. 

The paper is organized as follows. In Section [2] we introduce assumptions and 
notations. Section [3] includes three subsections concerning existence, uniqueness and con- 
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vergence. In Section [4] we present a method for the solution of the discrete system and 
some numerical examples illustrating the theory. 


2 Notations and preliminary definitions 


A network is a couple (V,£) given by a finite collection of vertices V := {uj}j g / and a 
finite collection £ := {ej}jeJ °f continuous non self-intersecting arcs whose endpoints 
belong to V. We assume that each arc ej G £ is parametrized by a smooth function 
TTj : [0. lj] —> M n , lj > 0. For a function u : T —> M we denote by Uj : [0, lj] —> M the 
restriction of u to ej, i.e. u(x) = Uj(y) for x G ej, y = nj 1 (x). Given G V, we denote 
by dju{vi) the oriented derivative at Vi along the arc ej defined by 


lim ( Uj(t ) — Uj(0))/t, if Vi = 7Tj(0); 

dju(v t ) ^ li m (uj(lj — t) — Uj(lj))/t, if Vi = nj(lj). 
t-¥ o+ 


Given a discretization step h = {hj}j<zj, we consider an uniform partition yj^ = khj, 
k = 0of the interval [0, l 3 \ which parameterizes the edge ej (we assume that 
Nj 1 = lj/hj is an integer). We obtain a spatial grid on T by setting 

Qh = {xj,k = 7T j(yj,k), j G J, k = 0,..., Nj}. (2.1) 

We define 


Inc+ = {j G Incj : Vi = vrj(O)}, Inc- = {j G Inc* : Vi = 7 tj(Njhj)}, 


so that 

as shown in Figure [l] 


Incj = Inc) 1- U InCj , 



Figure 1: incident edges to the vertex vf Inc) 1 " = { j }, InCj = {k,l}. 

We set 

|/!|=ma x{hj}, « , ‘ = #(/) + y(iV' 1 -l). (2.2) 

jeinci jeJ 

i.e. N h is the total number of the points of Qh, having identified for each i G / the ^t(lncj) 
grid points corresponding to the same vertex vy. For a grid function U : ^ M we 

denote by Uj t k its value at the grid point Xj t k- 
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Definition 2.1 We say that a grid function U : Q k —> M is continuous at Vi if 

Uj/ = U k , m := Ui if Vi = 7 Tj{thj) = ir k (mh k ), j,k G Inci, i G {O,^}, and m G {0 ,N%}, 

i.e. the value of U at the vertex v l is independent of incident edge ej, j G Inci. We say 
that a a grid function is continuous if it is continuous at Vi, for each i G I. 

We introduce the finite difference operators 

{D+U)jk = ViMl-UjJ. ' 

[D h U] jyk = (( D+U) jtk ,(D+U ) j , k - 1 ) T , 

! n2m Uj k— 1 — + Uj k +1 

\ 1J h U )j,k = - 72 -■ 

In order to approximate the Hamiltonian Hj : [0, If x E -> R, j G J, we consider a numer¬ 
ical Hamiltonian gj : [0. If X M 2 —> M, (x,qi,q 2 ) —> gj (x,qi,q 2 ) satisfying the following 
assumptions: 

(Gl) monotonicity: gj is nonincreasing with respect to q\ and nondecreasing with re¬ 
spect to q 2 . 

(G2) consistency: gj (x, q, q) = Hj(x , q) Vx G [0, lf\, Vq G M. 

(G3) differentiability: gj is of class C 1 . 

(G4) superlinear growth : gj(x, q\, qf) > a((qf) 2 + [qf) 2 ) 1 ^ 2 — C for some a > 0, C G R, 
7 > 1 and qf denote the positive and negative part of q s , s = 1, 2 
(G5) convexity : for all x G ej, ( qi,q 2 ) >->• gj (x,qi,q 2 ) is convex. 

Numerical Hamiltonians fulfilling these requirements are provided by Lax-Friedrichs or 
Godunov type schemes, see {22]. As an example, suppose that the Hamiltonian H is of 
the form H(x,p) = ^(x, \p\) where \k is convex, increasing and superlinear with respect 
to its second argument. Then the Engquist-Osher Godunov scheme reads as 

9j{x,qi,q 2 ) = T(x, (min(gi,0) 2 + max(g 2 ,0)) 2 ) 

and the monotonicity, consistency and coercivity conditions are satisfied. 


Given U,W : Qh — > M, we define the scalar product 
Nf-i / 

(U,W) 2 = J2 E h jU 3 , k W 3 , k + J2 I £ Y U hoWj,o+ Y y Uj, N ?W jiN H 

ieJ k =l iel \jeinc+ jeincr 

We introduce the compact and convex set 

K-h = {( M 3 ,k) je j. 0 <k<Nf '■ M is continuous, M^ k > 0, (M, 1) 2 = 1}. 

The operator V[m](xj jk ) is approximated by (V k [M])j k where M is the piecewise 
constant function taking the value Mj >k in the interval | y — yj lk \ < hj /2, k = 1,, Nf — 1, 
j G J (at the vertices only the half interval contained in [0, lj\ is considered). In particular, 
if V is a local operator, i.e. V[m\(x) = F(m{x )), then we set (Vh[M])j jk = F(Mj jk ). We 
assume that 
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(VI) Vh is continuous and maps JCh on a bounded set of grid functions. 
(V2) Vh is monotone, i.e. 


(V h [M] - V h [M],M - M) 2 < 0 =► M = M. 

(V3) There exists C independent of h such that for all grid functions M G JCh 

||14[M]|| 0O :=max\{V h [M]) jtk \<C 

j,k 

I (V h [M])j,k - {V h [M])jj I < C\y jtk - Vj/ \ M = 0,..., N}, j G J. 


3 A finite difference scheme for the stationary MFG system 


In this section we introduce the approximation scheme for the system (1.1). For simplic¬ 
ity, we consider a network T without boundary; appropriate boundary condition can be 
inserted in the scheme in a straightforward way. At the internal grid points we consider 
the finite difference system 

+ g( Xj , k , [D h U] j)k ) + A = (V h [M})- k , k = 1,..., N? - 1, j G J 


Vj {DlM) jik + B h {U,M) jik = 0, 
. M G JCh , (U, 1)2 = 0, 


k = I, - ■ ■, A/A - 1, j G J 


(3-1) 

where U, M are grid functions and A G M. The transport operator B h is defined for 
j & J and k = 1 by 


B h (U,M) jjk = 


Mj,kj^( x j,k, [DhU]j,k)+ 

Mj >k +lg&( x j,k+l, [Dh.U]j }k + 1 ) — A ^j,kQ§z{ x j,ki [DhU]j : k) 


for k = 2,..., Njf — 2 by 


B h (U,M) jtk = 


Mj,k Qq t ( x j,ki \PhJJ]j,k) Mj, k -l Qg ± { x j,k- 1) [DhU]j,k- 1) 


9g 


T Mj,k+1 Qq 0 { x j,k+li \DhU\j,k+ 1) Mj,kQq 2 ( x j,ki \DhfJ\j,k) 


dg 


for k = Nj — 1 by 


B h (U,M)j k = 


[DhU]j,k) — [PhU]j,k-\) — 


9g 


Mj,kdq 2 ( X j,ki [DhU\jjf) 


We discuss now the transition conditions at the vertices, see (1.2). We discretize the 
Kirchhoff condition for the function u via a 1 st order approximation of the derivative and 
we impose the continuity at the vertices 


S h {U,V h [M}-A)i = 0, 
U continuous at Vi, 


i G I, 
i G I, 


(3.2) 
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where for grid functions U, V, the operator S h : V —> M is defined by 

S k (U,V)i= £ h( D+ C)j,o + h^o] - •£ (3.3) 

ielncf jeinc- 

To discretize the transition condition for m we consider a 1 st order approximation of the 
derivative (the continuity of M at the vertices is included in the definition of Kh) 

T h (M,U)i = 0 i€l, (3.4) 

where for grid functions U, M, the operator T : V —> M is defined by 

T h (M,U)i= [D h U] jt i)] 

j£lnc+ 

da (3-5) 

— [PhU]j' N h_-^)\ = 0. 

jiGlnc” 


Remark 3.1 For the discretization of the differential equations in (1.1) inside the edge 


we follow the same approach in mw and we refer to these papers for motivations and 
explanations. We just recall that the transport operator B h comes from the 
discretization of the quantity 


/ mHpfx, du)dw dx 
Je-j 

for a test function w, which is connected with the weak formulation of the 
Fokker-Planck equation on the network. 


For the the approximation of the transition conditions in (1.2), we use 
a standard 1 st order discretization of the normal derivative of u and m with 
the sign depending if the vertex corresponds to either the initial point or the 
terminal one in the parametrization of the edge. The flux term in the Kirchhojf 
condition for m is approximated in a upwind fashion depending always on the 
orientation of the edge. Finally the additional term — A) in (3.2), 

which vanishes for h —» 0, is necessary to obtain the identity (3.22) which plays 
a key role in the uniqueness and convergence results. 

Note that at a vertex v t , we have respectively jf(Ina) values Uj and ff(Ina) 
values Mj, corresponding to the restrictions of these functions to the incident 
edges ej, j € /nc,. Since (3.2) and (3.4) gives jf(lnci) linear conditions, the 
value of U and M at v\ is univocally determined. 


Summarizing the approximation scheme for the stationary problem (1.1) is given by the 


(3.1)-(3.5). In the next subsections we study existence, uniqueness and convergence of the 


scheme. 


3.1 Existence 


We prove existence of a solution to (3.1)-(3.5) by a fixed point argument. We 


preliminarily need to prove existence, uniqueness and regularity for the first 
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equation in (3.1) with transition conditions (3.2), see Lemma 3.4 This result 
is obtained, as in the continuous case, by approximating the limit ergodic 
problem (3.15) with the sequence of problems (3.6), which contains a zero order 
term pU p , and passing to the limit for p — > 0. For this we need to estimate, 
uniformly in p, the discrete gradient of U p (see Lemma 3.14). 


Lemma 3.1 Let V : Qh —> M be a continuous grid function and assume that g satisfies 
(G1)-(G3). For p > 0, there is a unique solution to the problem 


~ u i( D h UP )i,k + 9( x j,k, [D h U p ] j>k ) + P U p . k = V j)k , k = 1,, Nh - 1, j € J 
S h (U p , V — pU p )i = 0, i€J; 


U p continuous at m, 


i E I. 


Proof To prove the existence we show that the map T : ~fil Nh —> M iV " defined by 
\ {fi.j( D h U )j,k ~ 9( x j,ki [D h U] jik ) + Vj,k ), j E J, k = 1,..., IVj 1 - 1; 


(3.6) 


pN h 


Hu) = 


±-S k (U,V)i, i e I; 


(where h Vi as in ( |2.2| )) admits a fixed point. 

Set r = (maxj^ \H(xj tk , 0)| + ||F||oo)//0- By the regularity of g the map T is contin¬ 
uous from B r = {U E M. Nh : ||C/||oo < r} to R Nh . Assume that U E dB r , hence 
max Jg j fc=0 N h\U j)k \ = r. Consider first the case Uj_ k = r for some j E J, k E 

{1 — 1}. Since (D k U)j jk < 0, D + Uj :k < 0 and D + Uj k _\ > 0, by the mono¬ 

tonicity and the consistency of g we get 


9{ x j tk ,[D/ l U]j^ k ) < H{xj tk , 0 ) 


and therefore 


HU)j, k ^ “ i~ H (xj, k , o) + Vj t k) < r. 


Hence F(U) jik < U jjk and F(U) j)k pU j)k if p > 1. 

Now assume that there exists i E / such that Ui = r for some i E / ([/* is the common 
value of Uj k at vf) then ( D + U)jp < 0 if Vi = ttj( 0), (D + U), N h_i > 0 if Vi = 7r ? (A r ( l ) and 

' J' j J 

therefore 

2 / hj „, v—v hi 


HU)* < 


ph Vi 


EE 


jelncf 




< r. 


jGlnc, 


Hence T(U)i < Ui and F(U)i pUi if p > 1. Arguing in a similar way if either Uj ik = —r 
or Ui = — r, we have that F(U) pU for all p > 1 and 1/ E dH r . Hence by the Leray- 
Schauder fixed point theorem there exists U p E B r such that F(JJ P ) = U p and therefore 
a solution of (3.6). We also have the estimate 


Halloo < -(max \H(xj k , ( 

P 3,k 


+ 


(3.7) 


We prove uniqueness of the solution to (3.6). Let U 1 , U 2 be two solutions of (3.6) and 
assume by contradiction that msocj k (Uj k — U 2 k ) = 5 > 0. Consider first the case that 
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there exists j £ J, k € { 1, ..., IVj 1 — 1} such that Uj^ — U?^ = 6. Subtracting the equations 
satisfied by U 1 and U 2 , we get 

-VjiDliU 1 - U 2 )) u r + g(x u , [.D h u\ k ) - g( X] - k , [iD h U \ k ) + p{U l - U\ k = 0. 
Since (j, k) is a maximum point for U 1 — U 2 , by the monotonicity of g, we get 

p5 = p(XJ l - U\ k < 0 

and therefore a contradiction. If there exists i £ I such that — U 2 = 6, then subtracting 
the transition conditions satisfied by U 1 and U 2 , we get 

0= ^ {v^D+iU 1 - U 2 ))^- J ^{p(U l - U 2 ) 3fi ) 

i6lnc+ 

- E (’'iV+w 1 - - %/< u ‘- u 2 )^) <-j E 

jelnc“ jelnci 


and therefore also in this case a contradiction. We conclude that U 1 < U 2 and we prove 
in a similar way that U 2 < U 1 . □ 

In the next lemma, we get an a priori bound for the gradient of the solution to the 
discrete Hamilton-Jacobi-Bellman equation by assuming that the function is bounded. It 
is important for the analysis of the convergence of the scheme that all the bounds are 
uniform in h. 

Lemma 3.2 Let V : Gh ~t ^ be a continuous grid, function and assume that g satisfies 
(Gl)-(GA). Let U h be a solution of the problem 

~ l 'j(D k U)j ! k + g(xj,k, [D h U]j,k) = Vj,k, k = 1,., N l j - 1, J e J 
< S h (U,V)i = 0 *€/; (3-8) 

U continuous at v^, i £ I, 

and assume that 

\\u h Hoc < Co (3.9) 

with Cq independent of h. Then 

pfc^lloo := max max \(D + U h )j k \ < C 
jeJ k=o,...,Nf-i 

where C depends on Co, [IV’Hoo, but not on h. 

Proof We first prove that D + U h is bounded at the vertices. Assume by contradiction 
that for some i £ I 

max < max \(D + U h )j : o\, nrax \(D + U h ) • N h_fi > —>• +oo for \h\ —> 0. 

I jelnct ’ je IncF ’ 3 \ 



Because of the transition condition in (3.8), it is not restrictive to assume that, up to a 
subsequence, 


max < max {(D + U h )j 0 }, max { — {D + U h ), N h_ 1 } > — > +00 for \h\ —> 0. 

[jeincf ’ jeinc- ■' J 

Hence we assume that there exists j £ Inc) 1 " such that D + Uj L 0 — > +00 for hj — > 0 (we 
proceed in a similar way if there exists j £ Inc)~ such that —D + U h N>} _ —> +00). 

Let ho be such that for hj < ho 

d + u^>-{c + \\v\u + ac ^ 

a Ij 

where C as in (G4), Co as in (3.9). Since (D k U h )j t i = [D + U^ l — D + Uj 0 )/h we have 


h 


- D+U h = TT D+U jfl + 9(x jt u [D h U\ 1 ) - V jtk > 
i n j 


^ D+U j ,o + <*\D + U^V - C - „. I|UU _ 


00 > ~r~^ + Uj t o 


and therefore D + Uj 1 > D + U ^ 0 . Iterating the previous inequality, we get 
D + U} )k+1 >D+U} >k for fc = 0,..., Nj — 1. 


For L < N-' — 1, we have 


U li = Ko + h 3 D+U^ 

Uf, 2 = Eft + = l/jjo + hj(D+u£ 0 + £+tft), 

EfcJ D + U* k - 


If Lhj > lj/2, by (3.11) we get 


Uh > Uh + LhjD+Uh > Eft + Lh,^ > Co 


J j,L ^ u j ,0 


j ,0 ^ u j ,0 


(3.10) 


(3.11) 


and therefore a contradiction to (3.9). 

We show that D + U h is bounded also inside T. Assume by contradiction that there 
exists j £ J, kh £ {1,..., — 2} such that, up to a subsequence, 


\ D+U l 


+00 


for h 0. 


(3.12) 


By compactness, Xj jkh —> xq £ for h 0. We set yo = A^o) £ [0, Zj] and we first 
consider the case yo £ (0, Zj). If D + U^ kh —> +00 for hj —> 0, let ho be such that for hj < ho 


Arguing as in (3.10), we have 


£> + tft fc+ i > D+U£ kh+1 _ 1 for / = 1,..., IVj 1 - A*. 
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For Lhj > (lj — yo)/2 we get 


L—l 


® + La+l = U* bh + h i Y, D + Ul kh+ , > -C 0 + Lh-^L- > C„ 


1=0 


lj - yo 


and therefore a contradiction to (3.9). 

We now consider the case yo G (0, lj) and D + U^ kh —> — oo for h —> 0. Let ho be such that 
for hj < ho 


owj^s-Lc+riu- — 

J ’ OL Vo 


We have 


T D+U i*. = (t d+ La-i + S(HU \D b u\ kh ) - v£ kh > 

Hj n 3 


+“i D+ ^kr - c - iLu > 


and iterating 
For L < kh — 1, we have 


0 + C* t „_, > D+Uj ,for / = 0,.... k h - 1. 


(3.13) 


La-i = La - 


Hence if Lhj > yo/2, by (3.13) we get 

L 


u ^ £ ^ - ^ D+E ^-' £ ^ - Lh > D+U ^ £ - c “ + ! c » 


/=! 


and therefore a contradiction to (3.9). 

In case yo = ' K J l (xo) is equal either to 0 or to lj and \D + U^ kh \ — > +oo, it is easy to 
adapt the previous arguments to obtain again a contradiction to (3.9). □ 


Lemma 3.3 Let U p be the solution of (3.6), then 

P + t/lco < C 2 

for a constant C 2 independent of p and h. 


(3.14) 


Proof Fix an arbitrary node Xj k G T and set W p = U p — U-- k - Adapting to the case 
of the networks the argument in (2[ Prop.2], it is possible to show that W p is bounded, 
uniformly in p. Since W p is a solution of (3.8) with V = V — pU p and by ( |3.7| ) V is 
bounded, uniformly in p and h, we can apply Lemma 3.2 to get a bound on ||LW p ||oo and 
therefore on ||ZDL/ p || uniform in p and h. □ 
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Lemma 3.4 Let V M be a continuous grid, function and assume that g satisfies 

(G1)-(G4). Then there exists a unique couple (U, A), where U : Q k —> M and A £ R, 
solution of the problem 

-"j(. D h U )j,k + g(x jjk , [D h U] jyk ) + A = V jtk , k = 1,..., Nh - l, j e J 
S h (U,V-A)i = 0, ie/; 

U continuous at Vi, i E I 

. (^1)2 = 0. 

Moreover 

|A| < Ci, II^CHoo < C 2 

for some constants C\, C 2 independent of h. 


(3.15) 


(3.16) 


Proof We prove existence by passing to the limit in the ergodic approximation (3.6). 

Halloo < Ci 


By (3.7) 


for any p > 0 where Ci is independent of p. By (3.7) and (3.14), up to a subsequence, 
U P — (U P , 1)2 converges to a function U : Qh ^ such that ( U , 1) 2 = 0 and pUj k converges 
to A E M (independent of (j, k)). Moreover the couple (U, A) satisfies (3.15) and the bounds 
in fl3l6l ). 

The uniqueness of the couple (U, A) can be proved by an argument similar to the one for 
the uniqueness of (3.6). □ 


Remark 3.2 Note that the dependence of the bounds in (3.16) on the function Vj yk is only 
by means of ||P||oo- This is crucial for the proof of the next theorem. 


Theorem 3.1 If g satisfies (G1)-(G4), V satisfies (VI), then the problem (3.1)-(3.5) 
has at least a solution (U,M, A). Moreover 

(3.17) 


|A|<Ci, ||t/||oo + ||^C||oo <c 2 

for some constants C±, C 2 independent of h. 

Proof We define a map <J> which associates to M 6 ICh the solution (U, A) of the 


problem (3.15) with V Jyk = (Vh[M])j yk . By Lemma 3.4 the map 4> is well defined. 

We show that is continuous. Let M s E ICh be such that M s —> M E ICh as s —> 00 , 
hence by (VI), Vh[M s ] Vh[M] as s —> 00 . Let (U S ,A S ) be the sequence of solutions 
of (3.15) with V = V h [M s }. By (3.16) the sequences A s and 


I U s 


are bounded and 


therefore, up to a subsequence, converge to A E M and, respectively, to a grid function U. 
It is immediate that (U, A) is a solution of (3.15) with V = Vh[M], By the uniqueness 
of the solution to (3.15), it follows that all the sequence (A s , U s ) converges to (A, U) and 


therefore the continuity of the map and the estimate (3.17). 


We define a map T which associates to Me ICh. the solution M of linear problem 


,Vj>-l 


pM jyk - u j (D 2 h M) jyk - B h (U , M) jyk = (j,Mj jk j E J, k = 1 

(•Mi * E iE ,„„+ ib(O + *k0 + + 

Ejelnc" i - i '/ 


i E I 
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where g > 0 and (U. A) = <h(M). We rewrite the previous problem as 


fiM + AM = nM 


(3.18) 


where A is N h x N h matrix. By the monotonicity and the regularity of g, for g sufficiently 
large the matrix gl + A is a non singular M-matrix and is therefore invertible. It follows 
that for any M G (3.18) admits a solution M and by M-matrix property M > 0 
since M > 0. We prove that (M, 1)2 = 1. First observe that if W, Z : Q}, —> M, then 


tv ?-1 


TV ?—2 


E HM d I w )m z m = -J2 E u i( D+w hA D+ z) jt 

j£j k =1 j£j k =1 


hk 


E E i z iAn+w)j »+E E 1 

jelncf 3 jelnc” 3 


and 


tv;- 1 


TV ?-1 


EE B h (U,W) jtk Z jjk = - EE bFj-fc • V q g{xj, k , 

jeJ fc=i jeJ fc=i 

-E( E ^M u z,,A m ,\ D h u] U )} 

*e/ jeinc+ J 

1 


ieinc^ 


If W = M, Z = 1, by ( |3.19[ )-( |3.20[ ) we get 

TV? —1 


El El (-AM) j,fc — El El j t . + Mj t \ q [DhU]j t i) 

jeJ k= 1 iS/ jelnc)- J 2 


hj L 


1 r 


E E E7 Z/ ?(-^ + -^ r )j,TVi l -l — Mj,N!?-l~flZ-( x j,N!}-l’ [^hU]j :N h_i) 


hi - 


jelnc i ^ 

Hence by the definition of .4 at the vertices we have 
tv;— 1 


<9<?i 


E E = ? E E - M t.» + E <* - M t. 


TV? 


jej fc=i 


\y_ !nc- 


Therefore 


tv;-i 


jelnCi 


tv;-! 


(3.19) 


(3.20) 


EE hj [gMj^ + {AM ) j^] — g EE hjMj tk 

j&J k=l j£j k =1 

which implies (M, 1)2 = (M, 1)2 = 1 and therefore M G /C^. 

Hence T maps K, k into /C^. From the boundedness and continuity of <1> and the regularity 
of g, \h is continuous. By the Brouwer’s fixed point theorem it follows that \h admits a 
fixed point M which is a solution of (3.1)-(3.5). □ 
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3.2 Uniqueness 

We first prove a fundamental identity which plays a crucial role in uniqueness and conver¬ 
gence of the scheme (compare with [3l (3.20)]. 


Lemma 3.5 Let A, B : Qh — 
and (U, M , A) a solution of 


be two grid functions, ( U,M,A ) a solution of (3.1)-(3.5) 


-Vj(D 2 h U) j)k + g{xj t k, [ D h U]j, k ) + A 

= {V h [M]). k + A j!k , 

k = 1,. 

Nb 

VjiDlM)^ + B h (U, M) jik = B jtk , 


k = 1,. 

Nb 

• • ’ E 


+ Ee/ncr 1 f A j,NA 

i G / 


T'‘(M, £?)i = + Y^anc- f 

i G I 


U continuous at v^, 


i G / 



. M G K h , (U, 1) 2 = 0. 


(3.21) 

Then 

n h {M, U, U) + n h (M, U, U) + (V h [M\ - V h [M],M - M) 2 = ( A, M - M) 2 + {B,U- U) 2 

(3.22) 

where 


Nb-i 


U h (M, U, U) = EE hjAlj}, \Dh,U\j,k) 

j&J k =1 

- [D h {U - U)\k,j ■ V q g{x jtk , [D h U] j}k ) 


Proof Let ( U , M, A) and (U, M , A) be as in the statement. Subtracting the equations 
for U and U, multiplying the resulting equation by hj(M — M)j }k and summing over j 6 J, 
k = 1,..., Nj 1 — 1, we get 


Nj~ 1 

E E hj{M-M) jjk 

j£j k =1 


VjDfafU U^j^k A \DhfJ\j,k) g(%j,ki [Dh,U\j,k)A 


(A - A) - iy h [M] - V h [M])j, h 


Nf-1 

EE hj(M - M) j>k A jjk . 

j&J k =1 


(3.23) 


Subtracting the equations for M and M, multiplying the resulting equation by h-j{U — U)jj ; 
and summing over j 6 J, k = 1 ,..., iVj ! — 1 , we get 




E E H u - A. 

j£j k =1 

Nf—1 

E E hj{U — U) jik B jt 

j£j k =1 


u D 2 h (M - M) jik + B h \U , M)j,fc - ^(U, M),-* 


(3.24) 


k* 
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N?-l 


We have the identity 

jyh-i 

-EE "j( U -U)j,k D l(M-M) j ,k = E E -Vj(M-M) jtk Dl(U-U) jik 

j£j k =1 jgj k =1 

+ E ^ [(M - M) il0 (D + (y - - ([/ - y) 3 ,o(^ + (M - M))^] (3.25) 

jeinc+ 3 

+ E i[(M - M) jtN H(D+(U - U))^ -(U- U) jtN H(D+(M - M))j^ 

je Inc- J 

and, respectively, 

JV^-l N h -2 

3 3 

E E B h (U,M) j}k (U - U)j, k = - E [ E M i* l D h( u - U)]j,k • [£> h C7] iifc ) 


jej fc=i 


j'eJ k=2 


- - U) jfi + Mj, N h_x( x j, N h- i > 


(3.26) 




In a similar way a corresponding equation for YljeJ Xk=L B h (U, M)j jk (U — U)j, k is also 
obtained. 

We now discuss the boundary terms in (3.25) and (3.26). By the transition conditions for 
U and U and the continuity of M at the vertices we have 


E £{M-M) jfl {D+(U-V)) jJ0 + E ^(M-M). Nh (D+(U-U)) 


jGlnct 


jglnc i 


_ h j 




3-1 


E 9 ( M - ^)j,o[(Vfc[M] - V h [M ]\0 - (A - A) - Aj-o] 


jelnct 


+ £ ^(M — M)^^[(I4[M] — Vh[M])- N h — (A — A) — Aj' N h]. 


jelnc^ 


By transition conditions for M and M we get 


(3.27) 


X] TT.(U ~ U)j,o(D + (M — M))j,o — E ^{U -U)^ h ( D+ {M - M)) jyNh _i 


i6lnc+ 


E 


{U-U)j S 


j elncT 


+ E 


jelnc. 


(tz-t/)^ 

hi 


hj 

j’GlnCj 

•i/,1 'N-O.l: [^^,1) - ID h U] jtl ) 


dg 


Q qi ( x j,N^-V i D hU] j>N h-i) 


o 1 

Mj,Nf-l~Q^( x j,N!?-l’ [DhUlj^h-x) + E 


j elnc+ 


+ E 0 B iM(U-U) 




jfG Inc- 


(3.28) 
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Replacing (3.25 )-( |3.28 ) in (3.24) and adding the resulting equation to (3.23) we finally 
get (recall that (M, 1)2 = (M, 1)2 = 1) 


K - 1 


^2 h i 

j£j k =1 


M. 


j,k 


[Dh,U]j t k) k) 

~[Dh{U - U)\ k j ■ V q g(xj tk , [. D h U]j tk ) 


+ Mj,k 


9ip^j,ki \BhU\j,k) 9(.3'j,ki \DhJJ\j,k) 


~[Dh{U - U)\ kJ ■ V q g{x jjk , [. D h U] j>k 
+ (V h [M\ - V h [M},M - M) 2 = (A,M- M) 2 + {B,U- U) 2 , 


which amounts to (3.22). 


□ 


Theorem 3.2 If g satisfies (G1)-(G5) and the operator V k is strictly monotone, i.e. 

(' V h [M] - V h [M],M - M) 2 < 0 => M = M (3.29) 


then the problem (|3.It) -(3.5) has at most one solution. 


Proof Let (17, M, A) and (17, M, A) be two solutions of (3.1)-(3.5). By (3.22) with 
A = B = 0 we get 

n h {M, U, U) + n h {M, U, U) + {V h [M\ - V h [M],M - M) 2 = 0. 

By the convexity of g and the monotonicity of V, we see that all the terms in the left 
hand side of the previous equality are positive and therefore must vanish. The strong 
monotonicity of V implies that M = M. Hence (17, A), (17, A) solve (3.15) with Vj :k = 

we get U = U and A = A. □ 


Vh[M]j ik = V h [M\ jjk and by Lemma 


3.4 


3.3 Convergence 


In this section we analyze the convergence of the scheme (3.1)-(3.5) in the reference case 

H(x,p ) = \p\ P + f{x) 


where /3 > 2 and / : T 


(3.30) 


is a continuous function. By [7], we know that in this case 


there exists a unique solution ( u,m , A) to (1.1) with u G C 2,a (r), m G C 2 (r), m > 0, and 
A G M. 

We consider a numerical Hamiltonian of the form 


g(x,p) = G(p 1 ,pj) + f(x) (3.31) 

where G(pi,p 2 ) = {p\ +p 2 ) 13 ^ 2 and pf denote the positive and negative part of p s , s = 1,2. 
We observe that g satisfies assumptions (G1)-(G5). We need an additional assumption: 

(V4) For any m G /C := {p G C 0,a (r) : f r gdx = 1}, M G K. h , denoted by Z h {M) the 
continuous piecewise linear reconstruction of M G lC h on T, then 

|| V[m\ - V h [M ]||oo < w(||m -X h (M)||oo) 

where u is a continuous, increasing function such that lim^ 0 + uj(t) = 0. 
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In the following we denote by o(l) a generic grid function whose maximum norm tends to 
0 as \h\ 0. Given a solution (u,m, A) of (1.1), we define a grid function u h by 


u j,k ■= T~. u j(y)dy, if j G J, k = 1,..., Nj - 1, 


dj J\y-yj,k\< h j/l 


■= E 

j'Elnc 


+ dj Jo<y-y j:0 <hj/2 


i(y)dy+ jr. [ 

. _ n .1 J 0 


,-eincr 


uj(y)dy, Hi el. 


(note that ( u h , 1)2 = 0). We define in a similar way the grid function m h E ICh and we 
also set X h := X. Observe that by (V4) 


lim || V[m\ - V k [m h ] ||oo = 0. 

h —>0 


(3.32) 


Hence by (3.32) and the consistency assumption (G2), ( u h ,m h , X h ) is a solution of 
' —v j {Df i u h ) j)k + g{x jtk , [D h u h ] jtk ) + X h = {V h [m h ]) jk + Aj >fc , fc = 1,..., N} - 1, j E J 
< + = fc = l I ... 1 JVf- 1 j€J 

. ^ G K-hi (u h , 1)2 = 0 
with the transition conditions 

' s>\ UK] - A'*)i = E, 6Imf ^ + e J 6 I .o- * € / 

■ ’ iel 

7 d 

u continuous at Vi, i E / 

where B h are two grid functions such that 


(3.33) 


(3.34) 


lim Halloo = 0, lim = 0. 

h—¥ 0 h—¥ 0 


(3.35) 


We need some preliminary lemmas. 


Lemma 3.6 Let (3 > 2, 

i) For all q, q E M 2 , 

g(x,q) - g(x,q)-V q g(x,q) ■ (q-q) > maxflp^ -2 , |p|^ _2 )|p - p| 2 (3.36) 

where p = {qf,qt), p = 

m) There exists a constant C such that for all q,q,r E M 2 and rj > 0 

l(V g s(s,9) - V q g(x,g)) -r\ < max(|p|^ 2 , |p|^ 2 ) ^|p - p| 2 + f?|r| 2 ^ , (3.37) 

where p= {qf,qt), P = (<?]"> 

For the proof of the previous lemma, we refer to [3j Lemma 3.2]. 
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Lemma 3.7 Let (U h ,M h ,A h ) be a solution of (3.1)-(3.5) and (u h ,m h ,X h ) a solution of 
(3.33)-(3.34) with m h > 5 > 0 for h sufficiently small. Then 


Nfy—i 


Jim E E h i l D+U kk-[D + Aj,k 


j£j k =1 


P 


= 0 


(3.38) 


Proof By the identity (3.22) with (U, M, A) = ( U h , M h , A h ) and (U, M, A) = (u h , m h , X h ) 
we get 

n h (M h , U h , u h ) + n h (m h ,u h , U h ) + [V h [M h ] - V h [m h ],M h - m h ) 2 

+ ( A h , M h - m h ) 2 + ( B h , U h - u h ) 2 = 0. 


By (3.17) and the regularity of u we get lim^o |(B , U h — u h ) 2 \ = 0. By m h ,M h € K-h 


and the Cauchy-Schwarz also get lining | ( A h , M h — m h ) 2 \ = 0. Hence by (3.36) we obtain 

(3.39) 


.Ad 1 


Sfe=i h j ml j,k max {MpE ' 3 2 }l p i!fc-Pj,fcl 2 = °(!)> 
j Ej£i hj Mh k max {| l^" 2 , |pj fc |^“ 2 } |Pj) fc - | 2 = o( 1), 


where P^ k = {{D + U h ) j k , (D + U h )f k _ 1 ) a nd p\ k = ((D+vf Y k , (D + u h )^ k _ 1 ). Since m h is 
strictly positive, by the first equation in (3.39) we get (3.38). □ 


Theorem 3.3 Let ( u,m,X ) be the unique solution of (1.1) and (U h , M h , A h ) the sequence 
of the solutions of the scheme (3.1)-(3.5). Then 

lim II U h - d|oo + || M h - m||oo + |A h - A| = 0. (3.40) 

Proof We set E h = M h — m h . Subtracting the equations satisfied by M h and m h and 
multiplying the resulting equations for we get 




E E [-MDlE% k + B h (U h ,M h ) jik -B h (U h ,m h )j, k 

jeJ k=l 


h A = 


Nf-l 


(3.41) 


E E [*V> m ' h kk - B h (U h ,m\ k + Bl k ] hjEj k . 

jaJ k =l 


By the transition conditions for M h and m h (recall that M h and, m h are continuous at 
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the vertices) 


E[ E E 'ifiH D+Eh )i,o+ E l^ 2 ^D h u%, 1 } 


*e/ j e in C + 


E E U P^W-i)] 


j’Elnc- 




<9gi 


*e/ jeinc+ 


k <9<?2 




E E j,N^\ m j,N^-l(g qi ^ X j,N^-l^[ D hU ]j :N h_ i) ^ (^AT^-i! [Dh.U ]j^ N h_ 1 )) + ^ B^ N h\ 


j'elnc. 


3<?i 


2 

(3.42) 


Arguing as in (3.19)-(3.20), we have 


N^-l JVj 1 -2 

E E D lE h ) j , k Ei k = -Y J E ^ip + ^ui 2 -E[ E ^iP + ^koi ; 

jeJ k= l jeJ fc=i *g/ jginc4 




+ 


E[- E ^.(^4.+ E g E E?-i(- D+E '‘t.»‘-. 


jiElnc^” 


j^Inc. 


+ e 

jglnc” 

Moreover 

JV?*-1 «t*-i 

EE [ B h (U h , M\ k - B h (U\ m\ k ] El k = -EE E\ k [D h E% k • V q g(x jtk , [D h U\ k ) 


j£j k =1 jdJ fc=l 

51 e *«**mww- E iw*A**^ 


jglnc. 


_ hj 


and 


^-i 

E E [B\U h ,m h ) j , k -B h (u h ,m h ) j>k 

j£j k =1 


E 


’j,k 


Nfy—i 


EE ]j,fe ‘ (V qS{Xj tk , [D k U ]j,k) ^ qd( x j,ki [DhU }j t k)) 

j&J k =1 


*eJ j e inc + 


jglnc 


_ 


3s- 


<9# 




<9f/i 


Set 


Nj‘—2 

(D h E h ,D h E h ) 2 =E E k^w+E [ E i( D+E E.oi 2 + E k d+e '‘W-iI 

jgj fc=l ■ - 


ie/ 


ieinc/ 


ieinc; 
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Replacing the previous equalities in (3.41), using (3.42) and recalling the estimate (3.16) 
we get 


K~ i 


(D h E h , D h E h ) 2 < -C [ £ E h AA,k + h A,k \DhE\ k ■ {V q g(x^ k , [iD h u\ k ) 

j£j k =1 

-V q g[ Xj ^[D h U\ k )) + Y,( E h i E j, 0 B j ,0 + E 

jElnc^" j'Glnc” 


(3.43) 


with C independent of h. By (|3.35), we have 


EmEE ■ E[ E h i E j,o B j,o + E <«a)(^^) 2 - (3-44) 

jElnc4 j'Glnc” 


Set P fe = ((P+Pj; fc )-,(P+P^_ 1 )+), / = ((P+uJ ;fe )-,(P+uJ jfe _i)+). By fU7]), ( ^ 
for any rj > 0 


m> j,k \ D hE h ]j, k ■ (V q g(x jtk , [DhU h ]j,k) ~ ^q9(xj lk , \D h U% k )) 


< A max(| Pl k f-\ | V h hk \A ~\ P lk ~ Pj,k I 2 + >71^* 


C, 


(3.45) 




U (jj\\D h U h U ~ lD h u h Uf + v\D k E^ . 


Plugging the estimates (3.38), ( |3.44 ), (3.45) in (3.43), we finally get 


(E h ,E h ) 2 + ( D h E\ D h E n ) 2 = o( 1) for \h\ -»■ 0. 

Hence we get the convergence of M h to m in P 1 (r) and uniform. By the convergence of 
M h to m and (V4), we get lim^o ||Ph[Af h ] —I4[ m,l ]lloo = 0. Hence U h and u h are solution 
of (3.15) with A = A h , Vj ik = Vh[M h ]j k and respectively A = \ h , Vj )k = V k [M h ]j tk + o(l). 
By a comparison principle for (3. 15[) , we get lA^ 1 — A ft | < o(l) and therefore 


lim |A — A"! = 0. 

|/i|-K) 


(3.46) 


Let u h be the continuous piecewise linear reconstruction of U h on T. By (3.17), u 


z-.h 


0, up to a subsequence. By (3.38) and (3.46), u is a weak solution to 
(1.1). Therefore by the uniqueness of the solution to (1.1), we get the convergence of U h 

□ 


uniformly as 
ter 

to u in H 1 ( r) and uniform. 


4 Numerical implementation and experiments 

This section is devoted to the implementation and test of a numerical solver for the sta¬ 
tionary MFG system (3.1)-(3.5). In [2], the stationary MFG system on the torus is solved 
via the so called forward-forward long time approximation : for a given approximation step 
h, the approximate solution (Uh, Mh, Ah) is obtained as the limit of (Ujf, Mfi, Uf/n/At) 
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for n —> oo, where (£/T, MV) is computed via discretization of the corresponding evolutive 
MFG system, implicit or explicit in time, up to time T = nAt. 

Here we propose a new approach which allows to compute the solution of the sta¬ 
tionary MFG system directly, avoiding long time or small delta approximations. We collect 
all the unknowns (U,M, A) in a single vector X oflength 2N h + l (with N h given by (2.2)) 
and we recast the 2N h + 2 equations of the stationary MFG system as functions of X. 
Hence we get a nonlinear map T : M 2Ari+1 —y M 2Ar ’ +2 defined by 

T h 
j 

rh 


X(X) = 


-Vj{D 2 h U)j )k + g(xj,k, [DhU]j,k ) + a - (Vh[M]) jjk k = 1 ,..., Nf - 1 , j E J 


v i{D\M) jik + B h (U, M) jtk 
S h (U, V h [M] - A)i 
T h {M , XJ)i 
( M , 1) 2 -1 
(U, 1)2 


k = 1, 
i G I 
i € I 




and we look for X* E 


1+1 such that 

F(X*) = 0 E 


p 2N h +2 


(4.1) 


By Theorem 3.1 and Theorem 3.2 there exists a unique solution to (4.1), but the system 
is formally overdetermined, having 2N h + 2 equations in 2N h + 1 unknowns. This ter¬ 
minology normally applies to linear systems, but is commonly adopted also in 
the nonlinear case with a slight abuse of notation. Indeed, the solution is meant 
in the following nonlinear-least-squares sense: 

X* = arg min ^ 11 J'(X) |||. 

To solve the above optimization problem, we employ the Gauss-Newton method, that we 
briefly recall here for completeness. We first denote the residual function by 

r(X) = \ WHX)Wl = \HX) T HX) 

and we consider the standard Newton method for approximating a critical point of r: 

U r (X k )6x = - Xr{X k ), X k+1 =X k + 5 x , k > 0, 
where the gradient Vr and the Hessian T~L r are given by 


■d 2 T t 


Xr(X) = J T (X) T X(X), U r (X) = J T (X) T JAX)+ -^(X)Fi(X), 


with 


BF 

(Mx)) iJ = zy-(x), 


BXi 


d 2 X t 

d 2 x 


{X) 


i —1 


d 2 X; 


k,e 


dx k dx t 


(X). 


Since we expect the residuals Ti{X k ) to be small for X k close enough to X *, it is reasonable 
to neglect the second derivatives in Fi r , using the approximation FL r {X) ~ Jjr{X) T Jjr(X). 
This yields the Gauss-Newton method: 


Jx(X k ) T JAX k )5x = -MX k Y HX k ) , x k+1 =X« + 5x, k> 0, 


rk\T 


k+1 _ 
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where the Jacobian Jj is well defined assuming that the numerical Hamiltonian g is of 
class C 2 in the gradient variable and that the operator Vh is of class C 1 . Despite this 
method allows to employ only first order derivatives of T, it is still not efficient from a 
numerical point of view. Indeed, once Jjr at X k is computed, we also need to assemble 
the right hand side jJj 7 and the matrix jJjj-, typically squaring the condition number 
of the system. This can be avoided by simply realizing that the k -th iteration of the 
Gauss-Newton method is just the normal equation for the following linear-least-squares 
problem 

mmh\J T (X k )6 x +HX k )\\ 2 2 , (4.2) 

Ox Z 

which is in turn easily and efficiently solved by means of the QR factorization of Jj. 
Indeed, let m = 2 N h + 2, n = 2N h + 1 and suppose that Jj^(X k ) = QR, where Q is 
a m x m orthogonal matrix (i.e. Q 1 = Q T ) and R is a m x n matrix of the form 

Ri 


R = 


0 


, with Ri of size n x n and upper triangular. Writing Q = {Q\ Qo) with 


Q i of size m x n and Q 2 of size m x (m — n), we get 


\\MX k )6x + X(X k )\\l = II Q t (MX k )6x + HX k )) 111 = II Q t QR5x + Q T X(X k 


R\$x 

0 


+ 


QlF{X k ) 

QlF(X k ) 


= \\Ri5 x + QjX(X k )\\ 2 2 + \\QlT(X k 


which is finally minimized by getting rid of the first of the two latter terms, i.e. solving 
the square triangular n x n linear system R\5x = — Q r {J-{X k ) via back substitution. 

Summarizing, we propose the following simple algorithm for the stationary MFG 
system: 

Given a guess A' = A 0 ), a tolerance e > 0 and a dumping parameter 

0 < a < 1 , 

REPEAT 

• Assemble F(X) and Jt{X) 

• Solve the overdetermined linear system Jx(X)5x = —F{X) in the least- 


squares SENSE (4.2), USING THE QR FACTORIZATION OF Jx(X) 

• Update X X + adx 

UNTIL 11 6x 112 < £ 

The algorithm is implemented in C-language and employs the library SuiteSparseQR m, 
which is designed to efficiently compute the QR factorization and the least-square solution 
to very large and sparse linear systems. 

Some remarks are in order: 

1) We always initialize the method by setting U° = 0, A 0 = 0 and M° = 1/L, where 
L = r lj is the total length of the network. In general there is no guarantee that the 


algorithm computes a minimum of (4.2) with zero residual, i.e. a solution of the stationary 


MFG system. Nevertheless, in all the tests performed, our algorithm seems to converge 
to a zero residual minimum independently on the initial guess. 

2) As for the standard Newton method, it is known that also the Gauss-Newton method 
may not converge if the dumping parameter is set to a = 1. A fine tuning of a. can 
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be accomplished via some moderate time consuming line search technique, but for our 
purposes we simply checked that the fixed value a = 0.9 is sufficient in all the considered 
examples. 

3) We never impose the constraint M > 0 in the computation. Surprisingly, our un¬ 
constrained optimization algorithm converges to a solution of the MFG system with non 
negative mass. We extensively checked this feature, also in the case of negative or changing 
sign initial guesses. Even if the mass can be negative in some intermediate iterations of the 
Gauss-Newton method, we always end up with a non negative mass in all the considered 
examples. 

4) Our technique can be successfully applied also in the homogenization of Hamilton- 
Jacobi equations, e.g. for computing the effective Hamiltonian for some cell problems. 
Our preliminary tests using the nonlinear-least-squares approach are very promising, both 
in terms of accuracy and computational costs. 

The previous points, in particular the convergence of the method, are still under 
investigation and will be addressed in a future work (see 0 )- 


We now set up the data for the numerical experiments. We consider a simple network 
in the plane with 2 vertices and 3 edges of unit length, as in Figure [2jr. For computa¬ 
tional purposes the network is mapped in a topologically equivalent network, 
in which one vertex is located at the origin and the edges are delimited by the 
3rd roots of unity (Vj = (cos(27rj/3),sin(27rj/3)), for j = 0,1,2), as in Figure [2^. 
Note that the boundary vertices, i.e. the vertices with a single incident edge, 
are identified and correspond to a single vertex on the network. 




Figure 2: a network with 2 vertices and 3 edges (a) is mapped in an equivalent network 
with boundary vertices identified (b). 


We assume that the numerical Hamiltonian has the form (3.31), with /3 = 2 and f(x) is 
such that, for j = 0,1, 2 and x G ej, 


f(x) = fj(x ) := Sj (1 + cos (27r(f + 1/2))) , x = tvj , t G [0,1] , 


where Sj £ {0,1} is a switch parameter to activate/deactivate the corresponding cost on 
the edge ej. If not differently specified, we discretize each edge by Nj = 250 nodes, so 
that the resulting nonlinear system has dimension 1502 x 1501, and we choose a tolerance 
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e = 10“ 4 for the stopping criterion of the algorithm. We finally assume a uniform diffusion 
on the whole network, i.e. Uj = u for j = 0,1,2 and v > 0. 

Here we are mainly interested in the qualitative behavior of the computed so¬ 
lutions, and we postpone at the end of the section some experimental analysis 
on the performance of the algorithm. Nevertheless, we remark that in all the follow¬ 
ing tests, the proposed method converges in about 10 iterations and the computational 
time is of the order of few seconds, even for larger grids. 

All the tests were performed on a Lenovo Ultrabook XI Carbon, using 1 CPU 
Intel Quad-Core i5-4300U 1.90Ghz with 8 Gb Ram, running under the Linux 
Slackware 14.1 operating system. 

Test 1. We consider a local operator of the form Vh [ rn ] = m 2 , and we choose a diffusion 
coefficient v = 0.1. Figure [3] shows the results corresponding to the activation of the cost 
/ on three, two or one edge, namely for (so,si,S 2 ) = (1,1,1), (so,si,S 2 ) = (1,1,0) and 
(sq,si,S 2 ) = (1,0,0) respectively. 


min = 0.039 max = 0.778 min = 5x10 4 max = 1.017 min = 0.053 max = 1.328 



Figure 3: the case Vh [rn] = m 2 and v = 0.1, the cost / is active on (a) three edges, (b) 
two edges, (c) one edge. 

In the top panels we represent the mass M using a color-map in which the blue and the 
red correspond respectively to the minimum and maximum values. Moreover, we repre¬ 
sent the network as a fatten tube, whose cross sections have a size proportional to M at 
the corresponding points. In the bottom panels we represent the network (in black) and 
both the mass M (in blue) and the corresponding value function U (in red). Since V is 
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increasing, it penalizes concentration of the mass. The cost / has, if sj 7 ^ 0, a maximum in 
the center of the edge e r Hence, if v is not to small, the agents should be well distributed 
on the network with a maximum of m around the minima of the value functions u, i.e. in 
the center of the edges where the cost is active. In fact we observe this behavior in all the 
three examples. 

Test 2. We are interested in the behavior of the solution as v —> 0, hence we choose the 
same parameters of the previous test, but with v = 10 -4 . In this respect, our method 
seems very robust and we can reach very small values of v even for quite coarse grids. 
Figure [4] shows the corresponding results. In this case we see that the solution is not 
better than Lipschitz and the support of Du and m are disjoint, as in the Euclidean case 
(see [2, Test 2]). 


min = 0 max = 0.939 min = 0 


max = 1.129 min = 0 max = 1.915 



Figure 4: the case Vj,, [rn] = m 2 and v = 10 4 , the cost / is active on (a) three edges, (b) 
two edges, (c) one edge. 


Test 3. We set V[m] = 1 — -arctan(m), we consider both v = 0.1 and v = 10~ 3 and 
the cost / active on the whole network, i.e. (so,si,S 2 ) = (1,1,1). Figure [5] shows the 
corresponding results. Since V is decreasing, the agents want the share the same position 
and therefore tend to concentrate around the minima of the value function. Note that for 
v small, the regularizing effect of the diffusion is small and m is close to a sum of Dirac 
functions concentrated at the minima of u. In this case assumption (3.29) is not satisfied 
and uniqueness of the solution may fail. 
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Figure 5: the case Vy [m] = 1 — ^ arctan(m) with (a) v = 0.1 and (b) v = 10 3 , the cost 
/ is active on the whole network. 


Test 4. In this experiment, we show that the method can efficiently handle the computa¬ 
tion on much more complicated structures. To this end, we consider the periodic network 
shown in Figure [6^. It is a self-similar set, in which the length of each edge scales with 
a factor 1/2 when moving to adjacent edges. Starting from the longest edges, we stop at 
the second level of branching and we identify the extremal boundary vertices. Moreover, 
we choose the local operator Vy [rn] = m 2 , uniform diffusion coefficient v = 0.1 and the 
cost / as before, active on the whole network. In this example the players are distributed 
on all the edges with a scaling factor which depends on the length of the edge. 

Performance and convergence. Here we present some results showing the con¬ 
vergence and performance of the proposed method, both in terms of accuracy 
and computational times. We consider the same setting of Test 1, with the 
cost / active on all the three edges of the network, i.e. (so,si,S 2 ) = (1,1,1). 
Moreover, we choose the same number of discretization nodes for each edge, 
namely Nj = N for j = 0,1, 2 and a variable N, so that the space step is h = 1/N 
on the whole network. Note that, in the present case, the total number of de¬ 
grees of freedom ( dofs ) of the problem is much more than N. Indeed, we have 
N nodes for each of the three edges and for both U and M, that is dofs= 6 N. 
Since the exact solution is unknown for this problem, we assume as correct 
the solution computed for N = 2000, denoted by (U ex ,M ex , A ex ). Then we define 
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min = 0 


max = 0.103 



(a) 




(b) (c) 


Figure 6: solution on a self-similar network, (a),(b) the mass M, (c) the value function U. 


the error as 

E h = || U - U ex ||i + || M - M ex ||i + |A - A ea j, 

where the discrete 1 —norm, for a generic vector W with 31V components, is 
computed as ||Wj|i = hYlk~i\^k\ and the exact solution is projected on the 
corresponding grid via linear interpolation. Finally, we define the experimental 
order of convergence as Eoc(hi, h^) = log(£V ll /£V l2 )/log(/ii//i 2 ) and we set e = 10 -8 
for the stopping criterion of the algorithm. 

In Figure [7 ^l we show, for N = 1000, the behavior of the computed A as 
a function of the number of iterations. In this case A ex = —1.058687, whereas 
A = —1.058876 is obtained after 20 iterations with |A — A ex \ = 0.000189. 

Similarly, in Figure [7fa we plot the error Eh for different space steps h, 
ranging from 10 -2 to 10~ 3 . This shows an experimental convergence at least of 
order 1. 

Finally, in Table [l] we report all the results, including the error |A — A ex \ 
related only to the approximation of the ergodic constant, the Eoc computed 
for successive space steps, the number of iterations and the corresponding 
computational times. We clearly see that, even for quite coarse grids, we get 
a reasonable approximation of the solution with a very low time consumption. 
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Iterations h 

(a) (b) 


Figure 7: convergence, (a) A vs number of iterations for N = 1000, (b) error Eh vs space 
step h. 


N 

Dofs 

Error Eh 

Error |A — A. ex \ 

Iterations 

Eoc 

Cpu time 

100 

600 

0.01159 

0.003737 

7 

- 

0.13 

200 

1200 

0.00544 

0.001734 

7 

1.09 

0.37 

400 

2400 

0.00241 

0.000762 

17 

1.17 

4.09 

800 

4800 

0.00091 

0.000284 

16 

1.40 

19.57 

1000 

6000 

0.00059 

0.000189 

20 

1.94 

47.94 


Table 1: Performance of the proposed method. 
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